Second-order rotational effects on the r-modes of neutron stars 
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Techniques are developed here for evaluating the r-modes of rotating neutron stars through second 
order in the angular velocity of the star. Second-order corrections to the frequencies and eigenfunc- 
tions for these modes are evaluated for neutron star models. The second-order eigenfunctions for 
these modes are determined here by solving an unusual inhomogeneous hyperbolic boundary-value 
problem. The numerical techniques developed to solve this unusual problem are somewhat non- 
standard and may well be of interest beyond the particular application here. The bulk-viscosity 
coupling to the r-modes, which appears first at second order, is evaluated. The bulk-viscosity 
timescales are found here to be longer than previous estimates for normal neutron stars, but shorter 
than previous estimates for strange stars. These new timescales do not substantially affect the cur- 
rent picture of the gravitational radiation driven instability of the r-modes either for neutron stars 
or for strange stars. 
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I. INTRODUCTION 

Recently the r-modes have been found to play an inter- 
esting and important role in the evolution of hot young 
rapidly rotating neutron stars. Andersson |l| first real- 
ized and Friedman and Morsink B confirmed more gen- 
erally that gravitational radiation tends to drive the r- 
modes unstable in all rotating stars. Lindblom, Owen, 
and Morsink B then showed that the coupling of grav- 
itational radiation to the r-modes is sufficiently strong 
to overcome internal fluid dissipation effects and so drive 
these modes unstable in hot young neutron stars. This 
result has been verified by Andersson, Kokkotas, and 
Schutz Q. This result seemed somewhat surprising at 
first because the dominant coupling of gravitational ra- 
diation to the r-modes is through the current multipoles 
rather than the more familiar and usually dominant mass 
multipoles. But it is now generally accepted that gravita- 
tional radiation does drive unstable any hot young neu- 
tron star with angular velocity greater than about 5% 
of the maximum (the angular velocity where mass shed- 
ding occurs). This instability therefore provides a natu- 
ral explanation for the lack of observed very fast pulsars 
associated with young supernovae remnants. 

The r-mode instability is also interesting as a possible 
source of gravitational radiation. In the first few min- 
utes after the formation of a hot young rapidly rotat- 
ing neutron star in a supernova, gravitational radiation 
will increase the amplitude of the r-mode (with spher- 
ical harmonic index m. = 2) to levels where non-linear 
hydrodynamic effects become important in determining 
its subsequent evolution. While the non-linear evolution 
of these modes is not well understood as yet, Owen et 
al. PI have developed a simple non-linear evolution model 
to describe it approximately. This model predicts that 
within about one year the neutron star spins down (and 



cools down) to an angular velocity (and temperature) 
low enough that the instability is again suppressed by 
internal fluid dissipation. All of the excess angular mo- 
mentum of the neutron star is radiated away via gravita- 
tional radiation. Owen et al. H estimate the detectabil- 
ity of the gravitational waves emitted during this spin- 
down, and flnd that neutron stars spinning down in this 
manner may be detectable by the second-generation ( "en- 
hanced") LIGO interferometers out to the Virgo cluster. 
Bildsten |Q and Andersson, Kokkotas, and Stergioulas 0] 
have raised the possibility that the r-mode instability 
may also operate in older colder neutron stars spun up by 
accretion in low-mass x-ray binaries. The gravitational 
waves emitted by some of these systems (e.g. Sco X-1) 
may also be detectable by enhanced LIGO S. Thus, the 
r-modes of rapidly rotating neutron stars have become a 
topic of considerable interest in relativistic astrophysics. 
The purpose of this paper is to explore further the 
properties of the r-modes of rotating neutron stars. The 
initial analyses of the r-mode instability |p|-p[ were based 
on a small angular-velocity expansion for these modes 
developed originally by Papaloizou and Pringle p| . This 
expansion in powers of the angular velocity kept only 
the lowest-order terms in the expressions for the vari- 
ous quantities associated with the mode: the frequency, 
velocity perturbation, etc. This lowest-order expansion 
is sufflcient to explore many of the interesting phys- 
ical properties of these modes, including the gravita- 
tional radiation instability. However, some important 
physical quantities vanish at lowest order and hence a 
second-order analysis is needed jl^]. For example the 
coupling of the r-modes to bulk viscosity vanishes in 
the lowest-order expansion. Estimates of this impor- 
tant bulk-viscosity coupling to the r-modes have been 
given by Lindblom, Owen, and Morsink pj, Andersson, 
Kokkotas, and Schutz |§,0, and by Kokkotas, and Ster- 



gioulas ||l^ . But none of these is based on the fully self 
consistent second-order calculation needed to evaluate 
this coupling accurately. Since bulk viscosity is expected 
to be the dominant internal fluid dissipation mechanism 
in hot young neutron stars, it is important to extend 
the analysis so that this important physical effect can be 
evaluated properly. 

The dominant internal fluid dissipation mechanism in 
neutron stars colder than about lO^K is thought to be a 
superfluid effect called mutual friction |l^ caused by the 
scattering of electrons off the magnetic fields in the cores 
of vortices. Levin ||lj] has shown that the importance of 
the r-mode instability in low-mass x-ray binaries depends 
crucially on the details of the mutual friction damping of 
these modes. Unfortunately the mutual friction dissi- 
pation also vanishes at lowest order in a small angular 
velocity expansion of the superfluid r-modes. Thus in 
order to evaluate this effect properly, it is also necessary 
to determine the structure of the r-modes of superfluid 
neutron stars through second order in the angular veloc- 
ity. This provides another motivation then for developing 
the tools needed to evaluate the second-order rotational 
effects in the r-modes. 

In this paper we develop a new formalism for exploring 
the higher-order rotational effects in the r-modes. Our 
analysis is based on the two-potential formalism |l5[| in 
which all physical properties of a mode of a rotating star 
are expressed in terms of two scalar potentials: a hy- 
drodynamic potential SU and the gravitational potential 
6^. We define a small angular velocity expansion for 
the r-modes in terms of these potentials, and derive the 
equations explicitly for the second-order terms. This ex- 
pansion provides a straightforward and relatively simple 
way to determine the second-order effects, such as the 
bulk viscosity coupling, that are of interest to us here. 
The equations that determine the second-order terms in 
the r-modes form an inhomogeneous hyperbolic bound- 
ary value problem that is not amenable to solution by 
standard numerical techniques. Therefore we have de- 
veloped new numerical techniques which could well have 
applications beyond the present problem. In particular 
these techniques will also be needed to solve the analo- 
gous superfluid pulsation equations that determine the 
effects of mutual friction on these modes. 

The timescales derived here for the bulk viscosity 
damping of the r-modes differ considerably from ear- 
lier estimates. We find the bulk-viscosity coupling to 
these modes to be weaker for normal neutron stars than 
any previous estimates. Consequently the gravitational 
radiation-driven instability is somewhat more effective 
at driving unstable the r-modes in hot young neutron 
stars than earlier estimates suggested. Although quanti- 
tatively different from earlier estimates, our new values 
for the bulk-viscosity damping time do not substantially 
alter the expected spindown scenario in hot young neu- 
tron stars. We re-evaluate the critical angular velocity 
curve (above which the r-mode instability sets in) and 
find no qualitative change from earlier estimates. Our 



new value for the minimum critical angular velocity is 
somewhat lower than earlier estimates: about 5% com- 
pared to about 8% of the maximum. In very hot young 
neutron stars there is the possibility that bulk viscosity 
could re-heat the neutron star (due in part to non-linear 
effects in the bulk viscosity) and so suppress the instabil- 
ity to some extent |16[ . This could result in a significant 
increase in the timescale required to spin down young 
neutron stars, and could therefore decrease significantly 
the detectability of the gravitational radiation emitted. 
Our new calculation of the bulk viscosity timescale in- 
dicates that reheating will not be a major factor in the 
evolution of young neutron stars. Our calculations also 
show that the bulk-viscosity coupling in strange stars is 
somewhat stronger than the initial estimates by Mad- 
sen [nTJ. We find that bulk viscosity completely sup- 
presses the r-mode instability in strange stars hotter than 
T ^ 5 X lO^K, in good qualitative agreement with Mad- 
sen. 

In Sec. H we review the structure of equilibrium stellar 
models through second order in the angular velocity of 
the star. In Sec. [II we review the two-potential formal- 



ism for describing the modes of rotating stars, and derive 
the small angular velocity expansion of these equations 
through second order. In Sec. IV we focus our attention 
on the "classical" r-modes, the modes found previously to 
be subject to the gravitational radiation driven instabil- 
ity. We obtain analytical expressions for the second-order 
corrections to the frequencies of these modes, and present 
numerical results for polytropes and for more realistic 
neutron star models. In Sec. M we develop the numerical 
techniques needed to find the second-order eigenfunctions 
for the r-modes; we use those techniques to find those 
eigenfunctions; and we present the results graphically. In 
Sec. VI we use our new second-order expressions for the 



r-modes to compute the effects of bulk viscosity on the 
evolution of these modes. In the Appendix we discuss the 
convergence of the numerical relaxation technique used 
in Sec. |V| to solve the unusual hyperbolic boundary value 
problem for the second-order eigenfunctions. 



II. SLOWLY ROTATING STELLAR MODELS 

Our analysis of the r-modes of rotating stars is based 
on expanding the equations as power series in the angular 
velocity fl of the star. The first step therefore in obtain- 
ing these equations is to find the structures of equilibrium 
stellar models in a similar power series expansion. This 
section describes how to solve the equilibrium structure 
equations for uniformly rotating barotropic stars in such 
a slow rotation expansion. The solutions will be obtained 
here up to and including the terms of order fl^. 

Let h{p) denote the thermodynamic enthalpy of the 
barotropic fluid: 
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(2.1) 



where p is the pressure and p is the density of the fluid. 
This definition can always be inverted to determine p{h). 
The barotropic equation of state, p = p{p), then deter- 
mines p{h) = p[p(h)]. The equations which determine 
the family of stationary, axisynimetric uniformly rotat- 
ing barotropic stellar models are Euler's equation, which 
for this case has the simple form 



= Va[h-'^rHl-p^)n'-<p], 
and the gravitational potential equation, 
VVa^" = -AirGp. 



(2.2) 



(2.3) 



In these expressions r and p ~ cos 6 arc the standard 
spherical coordinates, and $ is the gravitational poten- 
tial. 

We seek solutions to Eqs. ( ^.2[ ) and ( ^.3[ ) as power series 
in the angular velocity Q. To that end, we define 



h{r,p) = ho{r) + h2{r,p) 
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(2.5) 
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$(r,/i)-$o(r)+$2(r,M)^7^+0(^'), (2.6) 

where po is the average density of the non-rotating star 
in the family. Using these ex pres sions then, the first two 
terms in the solution to Eq. ( |2.2| ) are given by 



Co = /io(r)-$o(r), 



(2.7) 



C2 = /i2(^ m) - \T^Gpor\l - p^) - $2(r,/i), (2.8) 

where Co and C2 are constants. The non-rotating model 
can be determined in the usual way by solving the grav- 
itational potential equation, 
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(2.9) 



together with Eq. ( p.?]) . The integration constant, Co, 
can be shown to be Cq — —GMq/Rq by evaluating 
Eq. (2.7) at the surface of the star. The constants Mq 
and Ro are the mass and radius of the non-rotating star. 
The second-order contributions to the stellar struc- 
ture are determined by solving the gravitational poten- 
tial Eq. (2.3) together with Eq. (2.8). The second-order 



density perturbation p2 is related to /12 by 



P2{r,p) 



ft.2(r,/i). 



(2.10) 



Thus using Eq. ( p.8D , the equation for the second-order 
gravitational potential can be written in the form 



V'^Va$2+47rG 



dp 



$2 = 



47rG H^ {C2 + i7rGpor2[l-P2(p)]}, (2.11) 



dh 



wher e ^2/ ^) = ^(3^^ ^1)- We note that the right side of 
Eq. ( p.ll| ) splits into a function depending only on r plus 
a function of r mult iplied by P2 (p). Since the operator on 
the left side of Eq. ( 2.11 ) acting on P2{p) gives a function 
of r multiplied by P2 (p) , it follows that the second-order 
gravitational potential $2 rnust have a similar splitting: 



$2(r,p) = $2o(?-) + $22(r)P2(M). 



(2.12) 



Thus the partial differential equation ( ^.11 ) for $2 re- 
duces to a pair of ordinary differential equations for the 
potentials $20 and $22: 
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Appropriate boundary conditions are needed to selec t 
the unique physically relevant solutions to Eqs. (2.13) 



and (2.14). In order to insure that the gravitational po- 
tential is non-singular at the center of the star, r = 0, 
we must require that <i>20 and $22 satisfy the following 
boundary conditions there: 
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V dr 



= $22(0) 



(2.15) 



r=0 



The potential $22 must also fall to zero as r ^ 00. We 
can insure this by requiring that $22 match smoothly at 
the surface of the star to a potential that in the exterior 
of the star is proportional to P2{p)/r^. It is sufficient 
therefore to require that $22 satisfy the condition 
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(2.16) 



An additional condition is also needed to fix <i>20. It 
is customary to consider families of rotating stars which 
have the same total mass. In this case the monopole part 
of the exterior gravitational potential is the same for all 
members of the family. To ensure this, we must require 
that the potential $20 and its derivative vanish on the 
surface of the star: 



= $20(i?0) 



rf$2C 

dr 



(2.17) 
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It might appear that Eqs. (2.15) and (2.17) now over 
constrain the potential $20 ■ This would be the case, 
except that the constant C2 that appears on the right 
side of Eq. ( 2.13) is still undet ermin ed. The boundary 
conditions Eqs. (2.15) through (2.17) are just sufhcient, 
however, to fix uniquely the potentials $20 a-nd $22 to- 
gethe r wit h the integration constant C2 as solutions to 
Eqs. ( 2.13 ) and ( 2.14 ). We also note that these boundary 
conditions insure that 



r'^P2{r, i^jdrdfi. 



(2.18) 



In summary then, the thermodynamic functions h{r, fi) 
and p{r, f i) in slowl y ro tating barotropic stars are given 
by Eqs. (2.4) and (|2.5D, where p2{r, n) and h2{r, p,) are 
given by 



P2ir,p-) = 



^)/^('^'^^ 



|)jC2+<l>20(0 + i^GpV 

+ [$22(r)-i^Gpor']P2(A*)}- (2.19) 



These expressions for h{r, p) and p{r, p) depend only on 
the structures of the non-rotating star through the func- 
tions ho{r), Poj r), a nd {dp /dh) o, the potentials $20 and 
<I>22 from Eqs. (2.12) and (ETJ), and the constant C2. 

It is also instructive to work out an expression for the 
surface r = R{p, fi) of the rotating star. This surface oc- 
curs where the thermodynamic potential h[R{p, ^),p\ = 
0. Solving this equation, we find 



R{p,n) = Ro + R2{p) 
where R2{p) is given by 
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(2.20) 



i?2(/^)=i?20+i?22P2(M) 

-\C2 + \-nGpoR, 



AttGpqRq 
+ [^22iRo)-^7rGpoRl]P2{p)}. (2.21) 

We have developed a computer code that solves these 
equations numerically for stars with an arbitrary equa- 
tion of state. We have tested this code against analytical 
expressions which can be obtained for a polytropic neu- 
tron star equation of state, p = Kp^, with K chosen so 
that a I.4M0 model has a radius of Rq = 12.533km. We 
find that the constants that determine the slowly rotat- 
ing model for this polytropic case have the values C2 = 
.09802Co, R20 = .15198i?o, and R22 = -.37995i?o- Our 
numerical results agree with the analytical to floating- 
point precision. 



III. THE PULSATION EQUATIONS 

The modes of any uniformly rotating barotropic stel- 
lar model can be described completely in terms of two 
scalar potentials SU and 5$ 1l5| . The potential S^ is the 
Newtonian gravitational potential, while dU determines 
the hydrodynamic perturbation of the star: 



<5C/=^-<5$, 



(3.1) 



where Sp is the Eulerian pressure perturbation, and p is 
the unperturbed density of the equilibrium stellar model. 
We assume here that the time dependence of the mode is 
giujt and that its azimuthal angular dependence is e""*^, 
where lu is the frequency of the mode and m is an integer. 
The velocity perturbation (Su° is determined in this case 

by 



Sv" = iQ'^'VbSU. 



(3.2) 



The tensor Q"'' depends on the frequency of the mode, 
and the angular velocity of the equilibrium star: 



Q-b = 



1 



(w + mn)2 - 41^2 



{UJ + TOf7)(5"'' - 



4f]2 



u! + mQ 
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(3.3) 



In Eq. (3.2) the unit vector z° points along the rotation 
axis of the equilibrium star, S'^'' is the Euclidean metric 
tensor (the identity matrix in Cartesian coordinates) , and 
u" is the velocity of the equilibrium stellar model. 

In general, the potentials SU and (5$ are solutions of 
the following system of equations ||l^ : 

WaipQ^'^bSU) = -{uj + mn)^{6U + 5$), (3.4) 

ah 



dh 



(3.5) 



subject to the appropriate boundary conditions at the 
surface of the star for SU and at infinity for (5$. In order 
to discuss these boundary conditions in more detail we 
let S denote a function that vanishes on the surface of the 
star, and which has been normalized so that its gradient, 
TT-a = VqS, is the outward directed unit normal vector 
there, n'^Ua — 1- The boundary condition on the function 
6U at the surface of the star, E = 0, is to require that 
the Lagrangian perturbation in the enthalpy h vanishes 
there. Ah — 0. This condition can be written in terms of 
the variables used here by noting that 



Ah = Sh + 



5v°- 
ikVL 



^ah, 



(3.6) 



where k is related to the frequency of the mode by 



kQ 



mfi. 



(3.7) 



Thus using Eqs. (3J^) and (3^) the boundary condition 
can be written in terms of SU and 5$ as 



= 



KflidU + (5$) + Q'^'VahVbSU 



(3.8) 
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The perturbed gravitational potential (5$ must vanish at 
infinity, limr^oo'5*& — 0. In addition 5$ and its first 
derivative must be continuous at the surface of the star. 
The problem of finding the modes of uniformly rotating 
barotropic stars is r educe d th erefore to finding the so- 
lutions to Eqs. ( |3.4[) and (3.5) subject to the boundary 
condition in Eq. ( |3.8| ). 

The equation for the hydrodynamic potential SU, 



Eq. (3.4), has a complicated dependence on the frequency 
of the mode and the a ngu lar velocity of the star through 
Q"'' as given in Eq. ( |3.3| ). In the analysis that follows 
it will be necessary to have those dependences displayed 
more explicitly. To that end , we re-write Eq. (pA) and 
the boundary condition Eq. ( |3.8| ) in the following equiv- 
alent forms: 
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2mK 
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-K'{K^-^)n'^{5U + 5^), (3.9) 



f«;2,5'^''-4z'^z'') VaWf,(5[/ + '^:^^vd''V ah 6U 



= 0. 



(3.10) 



Here we use the notation w for the cylindrical radial coor- 
dinate, VJ = ryjl — fi^ , and -nj" to denote the unit vector 
in the w direction. 

Ou r pu rpose now will be to derive solutions to 
Eqs. (3.5) and ( |3.9| ) as power series in the angular ve- 
locity of the star. To that end we define the expansions 
of the potentials 5U and 5^ as 



5U 



Rln^ 



SUo + 5U2 



n^ 



nGpo 



d<i> = Rfp^ 



S^o + S^2 



n^ 



ttGpq 
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(3.11) 



(3.12) 



The normalizations of SU and (5$ have been chosen to 
make the SUi and 6^i dimensionless under the assump- 
tion that the lowest order terms scale as ft^. Here we have 
limited our consideration to the generalized r- modes [|8| : 
modes which are dominated by rotational effects and 
whose frequencies vanish linearly therefore in the an gula r 
velocity of the star. In this case k (as defined in Eq. 3.7) 
is finite in the small angular velocity limit, and so we 
may expand 



f22 

TrGpo 



(3.13) 



Using these expansions, together with tho se fo r the struc- 
ture of the equilibrium star from Eqs. (2.4) and (^.5|), 
it is straightforward to write down order by order the 
equations for the mode. The lowest order terms in the 
expansions of Eqs. (|3.9|) and (B^) are the following. 
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Po(«o'5" -4z"z'')V6(5f/o 



dp 



2mKr 



tu'^VaPo SUo = 0, 
(3.14) 



V"Va5$o = -4^G(^j (5C/o + <5$o), (3.15) 

Similarly, the lowest order term in the expansion of the 
boundary condition is 



{KlS''''-4z''z^)S/ahoVbSUo 
2mKo a 



H -w''\/aho SUo 
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(3.16) 



r=flo 



Continuing on to second order, the equations for the 
potentials are 
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V°Va^$2 = -47rG ( -jj- ] {SU2 + <5$2) 



dh 



-^vrGl^) h2{SUo + S^o)- (3.18) 



The second-order boundary condition is somewhat more 
complicated; it must include two types of terms. The 
first type comes from the second-order terms in the ex- 
pansion of Eq. (3.S) in powers of the angular velocity. 



The second type comes from the fact that the boundary 
condition is to be imposed on the actual surface of the 
rotating star, not the surface r — Ro- This second type 
of term is the corre ction to the lowest-order boundary 
condition, Eq. ( 3.16| ), needed to impose it on the actual 
boundary of the star (to second order in the angular ve- 
locity). Hence the terms of the second type are propor- 
tional to i?2, the second-order change in the radius of the 
star from Eq. (2.21). The resulting boundary condition 



+Ko(ko - 4)7rGpo('5C/o + (5$o) 
+i?2r'=V 



2mKQ 



w 



■SUow^Vahc 



0. (3.19) 
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In summary then, Eqs. ( 3.17]) an d (3.18) together with 
the boundary condition Eg. ( |3.19 ) determine the second- 
order terms in the structure of any generahzed r-mode. 



IV. THE CLASSICAL 7?-MODES 

There exists a large class of modes in rotating 
barotropic stellar models whose properties are deter- 
mined primarily by the rotation of the star |lq,M . We 
refer to these as generalized r-modes. In this section 
we restrict our attention however to those modes which 
contribute primarily to the gravitational radiation driven 
instability. These "classical" r-modes (which were stud- 
ied first by Papaloizou and Pringle P|) are generated by 
hydrodynamic potentials of the form (see e.g. Lindblom 
and Ipser ^) 



SUo = a[ — 



m+1 



Prn+l{f^)e" 



(4.1) 



It is str aight forward to verify that this 6U0 is a solution 
to Eq. (3.14) if the eigenvalue kq has the value 



Kq — 



(4.2) 



This SUq and kq also satisfy the boundary condition 
Eq. ( 3.16| ) without further restriction at the boundary 
(and at every point within the star as well). The grav- 
itational potential (5$o must have the same angular de- 
pendence as SUq. Thus, 6^0 must (through a shght abuse 
of notation) have the form 



,5$o = a<5$o(r)P™,i(A*)e""'^. 



(4.3) 



The gravitational potential Eq. (3.15) reduces to an or- 
dinary differential equation then for 5$o {t) '■ 



Once SUq and S^q are known, it is straightforward 
to evaluate the perturbations in other thermodynamic 
quantities to this order. For example Spo — podho ~ 
po{5Ua + S^Q). And it is straightforward to evalua te th en 
the velocity perturbation to this order using Eq. (3.2). 

We next consider the second-order contributions to the 
r-modes. First, let us analyz e the second-order equa- 
tion for the potential 5U, Eq. ( 3.17 ). This equation con- 
tains two types of terms: those proportional to 6U2, and 
those that are not. We will consider those terms not 
proportional to 6U2 as source terms, and we evaluate 
them now. It is convenient to break these source terms 
into three groups. The first group is proportional to p2- 
These ter ms can be simplified by recalling that SUq sat- 
isfies Eq. ( 3.14 ) for any spherically symmet ric d ensity 
distribution. Then, using the fact from Eq. (2. IE) that 
P2{r,fJ.) = P2air) + P22(?') ^2(^)1 we find 
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w^VaP25Ua 



12m(m + 2) P22 
(m -I- 1)2 r^ 



SUo. (4.5) 



The second group of terms is proportional to K2- These 
terms have the following simplified form: 

V"(2KoK2PoVa(5C/o) + ^^^n7°VaPo'5C/o = 

2(m + 2)«2-^<5[/o. (4.6) 
r dr 

Thus, combining toge ther these terms with those on the 
right side of Eq. ( 3.17 ), we obtain the following expression 
for the equation that determines 5U2 for the classical r- 
modes: 



Va< Pa 



4(5"'' 



Az^z^ 



yb5U2 



Amm^-VaPo 



(m + 1)2 

I2m{m + 2) P22 

{m + 1)2 r 



(m + l)w 

SUo-2{m + 2)K2-^6Uo 
z r dr 



SU2 



(4.7) 



A similar reduction can also be made on the second- 
order boundary condition, Eq. (2.19). We collect similar 
terms together to obtain the following simplifications: 
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r dr 
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TABLE I. The second-order eigenvalues K2 of the clas- 
sical r-modes for stars with polytropic equations of state 



n 


m = 2 


m = 3 


m = 4 


m = 5 


m = 6 


0.0 


.57407 


.59766 


.54720 


.49074 


.44044 


0.5 


.41718 


.43861 


.40415 


.36406 


.32782 


1.0 


.29883 


.32054 


.29946 


.27250 


.24729 


1.5 


.21183 


.23426 


.22369 


.20693 


.19019 


2.0 


.14777 


.17084 


.16846 


.15961 


.14942 


2.5 


.10091 


.12426 


.12808 


.12532 


.12016 


3.0 


.06716 


.09024 


.09859 


.10039 


.09905 


3.5 


.04334 


.06556 


.07699 


.08210 


.08357 


4.0 


.02692 


.04773 


.06102 


.06839 


.07186 


4.5 


.01589 


.03487 


.04896 


.05768 


.06252 



The latte r follows from the fact that the expressi on i n 
Eq. (pTfl) is zero everywhere if SUq is given by Eq. (4.1). 
Combining these simphfied expressions together gives the 
following form for the boundary condition that constrains 
5U2: 



ab 



12m{m + 2) /122 



Va/ioVbdC/2 + — — 6U2 

[m + l)vj 



6Uo + 2{m + 2)K2-^dUo 
[m + 1)-^ r-^ r dr 



16m(m + 2) 

(to +1)4 



irGpoiSUo + S^o) 



0. 



r=Ro 



(4.10) 



We note that the operator on the left side of Eq. (4.7) 
which acts on SU2 is identical to t he op erator that acts 
on 6U0 from the lowest-or der Eq . (3.14). We also note 



that the right side of Eq. (4/7) is a function of r multi- 
plied by the angular function P™_^]^(/i)e""'''. These facts 
allow us to derive a simple formula for the second-order 
eigenvalue K2 in terms of known quantities. Multiply the 
left side of Eq. ( [4.7|) by SUq and integrate over the in- 
terior of the star. This integral vanishes becaus e thi s 
operator is symmetric and 5U^ also satisfies Eq. (B.14). 
This implies that the integral of 6Uq multiplied by the 
right side of Eq. (4.7) also vanishes. This integral gives 
the following expression for the eigenvalue K2 once the 
angular integrals are performed: 

pRo / . \ 2m+2 , 
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(m + 1)2 7„ ^22 ^ ^^ 
STrGpom C^\ 

m-\-l 
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dr 



ni-\-l 



r 
Ro 



S^oir) 



dr. (4.11) 
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FIG. 1. Angular velocity dependence of the frequencies of 
the classical m = 2 r-modes for I.4M0 stellar models based 
on seven realistic neutron star equations of state. The dashed 
curve is based on the lowest-order expression for lu/Q, = kq — 2, 
while the solid curves are based on the second-order expres- 
sion LO/Q — Ko — 2 + K2i^'^ /llGpO- 



We have evaluated Eq. ( 1.11 ) numerically to deter- 
mine K2 for a variety of equations of state. Table | 
presents the values of K2 for the classical r-modes with 
2 < TO < 6 of stars with polytropic equations of state. 
We also present in Figure [| a graph of the frequency 
uj/Q = K — 2 of the TO = 2 classical r-modes computed for 
I.4M0 stellar models based on seven realistic equations 
of state pfl]. The dashed line in Fig. n^ corresponds to 
the lowest-order approximation of the r-mode frequency 
Lu/Q — kq — 2, which is the same for any equation of state. 
The solid curves are based on the second-order formula 
uj/n — kq — 2 + K2^^/ttGpo. It is interesting to see in 
Fig. H that the higher-order terms make only small (up to 
about 12% at the highest angular velocities) corrections 
to the frequencies of these modes for stars with realistic 
equations of state. We also note that the general ten- 
dency of the frequency of these modes to be smaller than 
that predicted by the lowest-order expression is consis- 
tent with the results found by Lindblom and Ipser |jl^ for 
the Maclaurin spheroids. An analytical expression for K2 
can be obtained from Eq. ( [4.11 ) for the uniform density 
case by performing the indicated integrals analytically. 
The resulting expression is equivalent to Eq. (6.10) of 
Lindblom and Ipser |18|. 



V. NUMERICAL SOLUTIONS FOR SU2 

In this section we discuss the numerical solution of 
the equations that determine the second-order correc- 
tions to the eigenfunctions 5U2 and 6^2 of the classical 



r-modes. Once 6U2 is known, the solution of Eq. ( 3.1§| ) 
to determine (5$2 is straightforward. Thus our discussion 
will concentrate on the more difficult problem of solving 



Eq. (4.7) for 5U2- It will be convenient to introduce the 
notation 



D{5U2)^Va\po 



4jab 



(m + l)2 



-Az'^z" 



S/b5U2 



4mzu"-VaPQ 
(m + l)w 



SU2, (5.1) 



for the diff erential oper ator that appears on the left side 
of Eq. (4/7). Thus Eq. ( li!?!) can be written in the form: 



D{SU2 



F, 



where 



12m(m + 2) p22^rr 

^ = ~f — rr^ — T^^o 

[m + ly r^ 



2{m + 2)K2-'^5Uo 
r dr 



-'«'<=^»i!^(l^i (»"+^*")^ (^'^ 



The problem of solving Eq. ( |5.2| ) numerically is a some- 
what non-standard problem that is made difficult by two 
facts. First, the operator D has a non-trivial kernel: 
D{SUo) = 0. Many of the straightforward numerical 
techniques fail in this case. Second, the operator D is 
hyperbolic. There appears to be little previous work on 
solving hyperbolic boundary value problems of this type. 
We solve Eq. (^) here using a variation of the stan- 
dard relaxation method commonly used to solve elliptic 
partial differential equations pl| . To that end we intro- 
duce a fictitious time parameter A and convert Eq. ( |5.2| ) 
into an evolution equation 



dxSU2 = D{8U2) - F. 



(5.4) 



The idea is to impose as initial data for Eq. ( p.4D a guess 
for 5U2, and then to evolve these data (as a function of 
A) until a stationary {d\5U2 = 0) state is reached. If suc- 
cessful, the late time solution (limA^oo 8U2) to Eq. (5.4) 
will also be a solution to Eq. (5.2). 



We implement the relaxation method to solve Eq. (5.2) 
by using a discrete representation of the functions and 
differential operators. Let u]-^ denote the discrete repre- 
sentation of the function 5U2 evaluated at the fictitious 
time A„. The index i (and later j and k as well) takes 
on values from 1 to A^ where N is the dimension of the 
particular discretization used. Similarly the discrete rep- 
resentation of the differential operator D of Eq. (5.2) is 
denoted D-'i, and the representation of the right side of 
Eq. ([5.2|) is denoted F^ . Thus the discrete representation 
of Eq. (|5.2) is simply 



D^,u' = F^. 



(5.5) 



Summation is implie d for pairs of repeated indices (e.g. i 
on the left side of Eq. 5^). The algebraic Eq. (5_^) cannot 
be solved by the most straightforward direct numerical 
techniques because the operator D has a nontrivial kernel 



(mentioned above). Consequently the matrix D^i has no 
inverse. 

Thus we are lead to introduce the evolution Eq. (pA). 
We use the "implicit" form of the discrete representation 
of Eq. (U): 



Oh<+i ^ [ D=^ - -^P^ <+i = F^ - ^<. (5.6) 



In Eq. ( |5.6| ) I^ i is the N dimensional identity ma- 
trix, and AA is the relaxation timestep. Given u\ we 



(5.2) solve Eq. (M) for 



n-l-l 



by direct solution of the lin- 



ear algebraic equation. We use the band-diagonal lin- 
ear equation solver LINSIS from EISPACK to compute 
[O-^YjiF^ - </AA). We find that this can be com- 
puted stably and accurately for almost any value of the 

relaxation timestep AA. 

Unfortunately solving Eq. (|5.6|) iteratively does not 
yield the desired solution to Eq. (|5.5| ) in the limit of large 
n. Instead the solution grows exponentially, becoming in 
the limit of large n closer and closer to a non-trivial so- 
lution to the homogeneous equation, 5Uq. Fortunately, 
this malady is easily corrected. Let -u* denote the discrete 
representation of 5Uq] thus D^iU^ « since SUq is in the 
kernel of D. Also we let Ui denote the discrete represen- 
tation of the CO- vector associated with u'. In particular 
choose Ui so that u'ui is the discrete representation of 
the integral of |(5C/oP- Then the matrix 



P',.^P. 



-k - ' 



(5.7) 



is the discrete representation of the operator that projects 
functions into the subspace orthogonal to 6 U0. We use 
this projection in conjunction with Eq. (5.6) to define a 



modified relaxation scheme to determine 1 



'n+l 



^n+1 



p'^o-ykiF'--^ 



AA 



iteratively: 



(5.8) 



By applying the projection P^j after each relaxation step 
we insure that the exponentially growing kernel is re- 
moved from the solution. We find that the iteration 



scheme defined in Eq. (5.8) does converge quickly and 
stably to a solution of Eq. ( p3.5D . In the Appendix we dis- 
cuss the reason this numerical relaxation method works 
even in the case of the unusual hyperbolic boundary- 
value problem considered here. We show that conver- 
gence is guaranteed for sufficiently large values of the 
relaxation timestep AA, and that it also converges for 
either sign of AA. 

In order to implement this inversion scheme we need 
explicit discrete representations of these operators. We 
find it convenient to work in spherical coordinates r and 
/i = cos 9. In terms of these coordinates then, the differ- 
ential operator D has the form 



D{SU2 



4po 



(m + l)2 



d^6U2 1 - M d^SU2 2 9(5[/2 
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r drd^ r^ dji^ 
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r 9r r 

4 f dp\ dho 



d^ 



[l-(m+l)V 



2, ,2 



9^ 

dr 



{rn + 1)2 \dh J Q dr 

(m I i')2 A*(l-M ^)^'^^2 , _______ , -,x"'^2 



,(5C/2 
+ m(m + lj 



(5.9) 

A similar spherical representation is also needed for the 
boundary condition, Eq. (3.11), 



2, ,2 



l~(m + l)> 



dho d5U2 , , , 1^1'^^o .^r 

— ; t: h mim + 1) —dU2 

dr or r dr 



( , Tx2 n 2^'^<^hQd5U2 
-(m + 1) m(1-M )-- 



1/ I 1^2 1 dh^ h22 




(m + 1) 



SUo 
= 0. (5.10) 



r=_Ro 



The operators in Eqs. ( |5.9| ) and (5.1C) are transformed 
into the discrete matrix representation of the operator 
D'^j using the techniques discussed in Ipser and Lind- 
blom ||l5[]. In particular we use a grid of points {rj,iik) 
where the rj are equally spaced in the radial direction 
and the fik are the zeros of one of the odd-order Legen- 
dre polynomials J2^ . We use standard three-point differ- 
ence formulae for the derivatives in the radial direction, 
and the higher-order multi-point formulae for the angular 
derivatives described in Ipser and Lindblom |lj]. Using 
this discretization of the operator D, we find that the it- 
eration scheme described in Eq. (5.8) converges rapidly. 



We begin the iteration by setting Ug = and find that 
after about five steps with AA = — 10^ Rq / pc (where pc is 
the central density) the changes in u!^ from one iteration 
to the next become negligible. 

Since the eigenfunction dU2 is somewhat complicated 
we present several different representations of it graph- 
ically. Fig. g depicts the functions SU2{r, pk) with ^k 
located at the grid points used in our integration: the 
roots of Pig{pk) = in this case. We present in Fig. g 
another representation of this dU2 in which we graph the 
functions SU2{rk,p) for r^ = -^Ro- This graph gives a 
clearer picture of the angular structure of 6U2- Finally 
we give in Fig. ^ another representation of the function 
6U2 in which we decompose the angular structure of 6U2 
into spherical harmonics by defining the functions fk {r) : 



5U, 




FIG. 2. Functions SU2{r,iJ.k) for a range of values of the 
angular coordinate /ik. The numbers along the right vertical 
axis are the values of k. These range sequentially from fii = 
at the equator of the star, to /xio « 0.992 near the rotation 
axis. The equation of state is the polytrope discussed in the 
text. 
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FIG. 3. Functions SU2{rk, A*) for a range of values oirk/Ro- 
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FIG. 4. Functions fk{r) that determine the sphe rical har- 
monic decomposition of 5U2 as defined in Eq. (5.11). 



<Jt/2(r,/i) 
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fkir)P, 



m+2fe-l- 



(5.11) 



We find numerically that the fk{r) are negligibly small 
except for the smallest few values of k. In Fig. Q we graph 
the first three /fc(r). 

To measure the degree to which 6U2 satisfies the orig- 
inal differential equation, we define 



/ \D{dU2) - F\^r^drdn 



(5.12) 



We find that the value of e achieved by a given solution 
is approximately e w (A.S/Nr)'^ where Nr is the number 
of radial grid points used in the discretization p3| . It is 
instructive to compare e to the quantity 



eo 



J\D{SUo)\^r^drdfi 
/ \F\^r^drdn ■ 



(5.13) 



which measures the degree to which SUa is in the kernel 
of D. Since analytically D{5Uq) = 0, the deviation of 
eo from zero is a measure of the accuracy of our discrete 
representation of D. We find that eo is approximately 
eo « [l.T/NrY in our numerical solutions. This scaHng 
is what is expected from the truncation errors involved 
in the three-point difference formulae used to construct 
D]. Since e < eo in our numerical solutions, we see that 



SU2 is a good solution to Eq. (5.2) 



VI. BULK VISCOSITY TIMESCALES 



E=h 



{pSv''6v*^ + dU6p*)d^x. 



(6.3) 



Bulk viscosity causes the energy in a mode to decay 
(or grow) exponentially with time. We can evaluate the 
imaginary part of the frequency of a mode th at re sults 
from bulk-viscosity effects by combining Eqs. (6.2) and 
( |6.3[ ). The result, which defines the bulk- viscosity damp- 
ing time, Tg, is given by 



1 

TB 




(6.4) 



In order to evaluate 1/tb we need to have explicit ex- 
pressions for th e va rious terms that appear in the inte- 
grands of Eqs. ( |6.2| ) and (3.3). The energy E, for exam- 
ple, can be expressed as the integral 



E = ^(m + l)3(2m + ly.R^n^ 
2m 



Mr) I ^ 



2m+2 



dr + Oin"^), (6.5) 



by performing the angular integrals indicated in Eq. (6.3) 



In the dissipation integral, Eq. (3.2), an explicit ex- 
pression for the bulk viscosity C is needed. In standard 
neutron-star matter the dominant form of bulk viscos- 
ity is due to the emission of neutrinos via the modified 
URCA process 25 1. An approximate expression for this 
form of the bulk-viscosity coefficient is pfl] 



One of our primary interests in evaluating these modes 
through second-order in the angular velocity is our desire 
to obtain the lowest-order expression for the bulk viscos- 
ity damping of these modes. Bulk viscosity is driven by 
the expansion, da = VaJw", of the fluid perturbation 
which can be expressed in terms of the scalar perturba- 
tion quantities 5U and S^ as 



Q'"'\7ah\7b5U + Kn{5U + S<P) . (6.1 



J. i dp 

^^ = — 7777 
p an 



The first term on the right side of Eq. ( |6.l| ) vanishes at 
lowest order. Thus the lowest-order contribution to the 
expansion 5a is at order fi'^. And thus the second-order 
quantities 5U2, H2 etc. that we have evaluated in the 
preceding sections are needed to evaluate 5a. 

Bulk viscosity causes the energy associated with a per- 
turbation to be dissipated according to the formula. 




-- = - / C5a5a*d^x 



(6.2) 



where C is the bulk viscosity coefficient, and E is the 
energy of the perturbation as measured in the co-rotating 
frame of the fiuid. The energy E can be expressed as an 
integral of the fluid perturbations: 



C = 6.0 X 10 



-59 P 



K2f]2' 



(6.6) 



where all quantities are expressed in cgs units. For the 
case of the classical r-modes the expansion 5a that ap- 
pears in Eq. (6.2) can be expressed explicitly in terms of 
the potentials 5Uo, (5$0i and 5U2'- 
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It will be of some interest to evaluate the accuracy of 
some previously published approximations for the bulk- 
viscosity timescale. One of these Q is based on an ap- 
proximation 5a « 5s for the expansion of the mode, 
where 



10 



Ss= 



p 

2i 1 (dp 
m + 1 po V '^''■ 



{6Uq + 5^a)Rln'^ + o[n^). 

(6. 



We note that 5s i s ju st the last term in the expres- 
sion given in Eq. ( p-TJ ) for 5a. It is the only term in 
Eq. ( |6.7D that depends only on the lowest-order perturba- 
tion quantities: the others depend on higher-order correc- 
tions through 5U2, K2 or /i22- We define the approximate 
bulk-viscosity timescale r^ i n an alogy with Eq. (6.4) by 
replacing 5a with 5s in Eq. (3.2). 

The bulk-viscosity contribution to the imaginary part 
of the frequency, 1/tb, is proportional to il^. This follows 
from Eqs. ( |6.2| ) and (6^) because E scales as fP from 
C, as n~^ from Eq. (6.6), and 5a as ft^ from 

||2^. The bulk-viscosity damping time, I/tq, 

also scales with temperature as T^. Thus it is convenient 
to define tb'- the bulk viscosity timescale evaluated at 
rj2 = TrGpo and T = 10^ K, 



Eq. (3^ 
Eq. (3.7 
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(6.9) 



We have evaluated tb numerically for the m — 2 r-mode 
(the one most unstable to gravitational radiation) of a 
1.4Mq stellar model with the polytropic equation of state 
discussed in Sec. ||, and find tb — 2.01 x lO^^s [p8[ . 
For comparison, we have re-evaluated the approximate 
timescale Tg described above and find f^ — 7.04 x 
lO^s " 
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The most interesting application of these dissipative 
timescales is to use them to evaluate the stability of ro- 
tating neutron stars to the gravitational radiation driven 
instability in the r-modes 0. The imaginary part of 
the frequency of the r-modes, 1/r, includes contributions 
from gravitational radiation tqj^ in addition to shear fs 
and bulk tb viscosity effects. The general expression for 
t{Q,T), a function of the temperature T and angular 
velocity fl of the star, is given by 



1 



1 



T{n,T) TGR 



n^ 



1 /10«K 



TS 



nGp 

TB vio^k; 



T 
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(6.10) 



The bulk viscosity timescale tb — 2.01 x lO^^s has been 
evaluated in this paper, while the gravitational radia- 
tion timescale, tgr = —3.26s, and the shear viscosity 
timescale fs = 2.52 x IQ^s, were obtained by Lindblom, 
Owen and Morsink ||^ for the polytropic stellar model 
discussed in Sec. O. It is interesting to determine from 
this expression the critical angular velocity flc'. 
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FIG. 5. Critical angular velocities for the r-mode instabil- 
ity as a function of temperature. The dashed curve gives the 
critical angular velocities based on the approximate bulk vis- 
cosity damping time Tg. 



For stars at a given temperature, those with Q > flc are 
unstable to the gravitational radiation driven instability 
in the r-modes, while those rotating more slowly are sta- 
ble. Fig. ^ depicts fie for a range of temperatures relevant 
for hot young neutron stars. For stars cooler than about 
lO^K, superfiuid effects change the dissipation processes 
completely and the analysis presented here is no longer 
relevant. The dashed curve in Fig. @ represents the crit- 
ical angular velocities computed using the approximate 
bulk viscosity damping time fg = 7.05 x lO^s instead of 
the exact value. We see that even though tb ~ 29f^, 
the qualitative shape of the flc curve is not affected. The 
minimum of the flc{T) curve depicted in Fig. ^ occurs at 
minilc = .OSOl^TrGpo which is 4.51% of the approximate 
maximum angular velocity ^^irGpo- 

Stars composed of strange quark matter are subject to 
a different form of bulk viscosity caused by weak inter- 
actions that transform d quarks to s quarks. The bulk 
viscosity coefficient that results from this process is given 
approximately by pO(| 



C = 3.2x 10= 



pT' 



( 



V 



K^n^ VlOOMev/ 



(6.12) 



where all quantities (except nis , the mass of the s quark) 
are given in cgs units. We have used this form of the bulk 
viscosity to estimate the damping time of the r-modes in 
strange stars. We find that in strange stars this damping 
time scales as 



1 _ 1 / T Y 



n^ 



nGpaJ VlOOMev 



(6.13) 



Using the polytropic stellar model described in Sec. O, 
we have computed the bulk-viscosity damping time of 
the r-mode to be tb — 0.886s for strange stars. This 
value is smaller (by about a factor of 7) than that found 



11 



by Madsen |17|, who used a very rough estimate of the 
bulk-viscosity damping time [|l2| for the r-modes. We 
also note that our expression for tb does not scale with 
angular velocity in the same way as Madsen's. Neverthe- 
less, our calculations confirm Madsen's prediction that 
bulk viscosity completely suppresses the r-mode insta- 
bility in hot strange stars. Using our expression for tb 
we estimate that the r-mode instability is suppressed in 
all strange stars with T > 5 x lO^K. 
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da = it is easy to choose AA so that \xa\ < 1 for all a, 
e.g., by taking A A sufficiently large. Thus, the sequence 
Un+i converges to 



lim u" 



XaAXF°' 
1 + x„ 



pa 

d„ 



(A5) 



Thus the implicit relaxation scheme converges to the de- 
sired solution to Eq. (5.2). 

In con trast to the implicit relaxation method defined in 
Eq. (5^), the analogous explicit relaxation method does 



not converge at all for this problem. An analysis similar 
to that carried out above reveals that the criterion for 
convergence of the explicit scheme is that jAAda + ll < 1. 
Clearly this can only hold for operators D where the 
eigenvalues all have the same sign (as is the case when 
D is elliptic) and only when AA has the correct sign. 
Our limited experience with hyperbolic operators D is 
that their eigenvalues have both signs. Consequently it 
is not surprising that our attempts at explicit numerical 
relaxation fail in this case. 



APPENDIX: NUMERICAL RELAXATION 



The operator D defined in Eq. (5.1) is symmetric, in 
the sense that 



g*D{f)d\ 



rD{g)d\ 



(Al) 



for arbitrary functions / and g in any stellar model where 
Po = on the surface. Thus the discrete representation of 
the operator Dij will be a Hermitian matrix, and conse- 
quently D^ j will have a complete set of eigenvectors. Let 
e^ denote the eigenvector corresponding to eigenvalue da : 
D^jB^^ = dacl^. Since these eigenvectors form a complete 
set, we can express any vector as a linear combination of 
them. Thus we take F^ = J^a ^"^a^ "Ki = Ha ""n^L' etc. 
The numerical relaxation scheme indicated in Eq. (5.S) 



can be re-expressed therefore in the eigenvector basis as: 



u. 



_ AXF" - < 
"+i AAd„-l 



The role of the projection operator P^j is merely to re- 
move from Eq. ( |A2| ) the component corresponding to the 
zero eigenvalue. The recurrence relation Eq. (A2) can be 
solved analytically: 



.^+l=x„AAF"^(-x„)^ 



where 



k=0 



1 



AXda - 1 



(A4) 



This series converges as long as \xa\ < 1- Since the pro- 
jection operator has eliminated the one equation where 
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